home *** CD-ROM | disk | FTP | other *** search
/ Language/OS - Multiplatform Resource Library / LANGUAGE OS.iso / lisp / elk-2_0.lha / elk-2.0 / src / bignum.c < prev    next >
C/C++ Source or Header  |  1992-10-09  |  17KB  |  713 lines
  1. #include <math.h>
  2. #include <ctype.h>
  3.  
  4. #include "scheme.h"
  5.  
  6. Object Make_Uninitialized_Bignum (size) {
  7.     Object big;
  8.     
  9.     big = Alloc_Object ((sizeof (struct S_Bignum) - sizeof (gran_t)) +
  10.            (size * sizeof (gran_t)), T_Bignum, 0);
  11.     BIGNUM(big)->minusp = False;
  12.     BIGNUM(big)->size = size;
  13.     BIGNUM(big)->usize = 0;
  14.     return big;
  15. }
  16.  
  17. Object Copy_Bignum (x) Object x; {
  18.     Object big;
  19.     register size;
  20.     GC_Node;
  21.  
  22.     GC_Link (x);
  23.     big = Make_Uninitialized_Bignum (size = BIGNUM(x)->usize);
  24.     BIGNUM(big)->minusp = BIGNUM(x)->minusp;
  25.     BIGNUM(big)->usize = size;
  26.     bcopy ((char *)BIGNUM(x)->data, (char *)BIGNUM(big)->data,
  27.       size * sizeof (gran_t));
  28.     GC_Unlink;
  29.     return big;
  30. }
  31.  
  32. Object Copy_S_Bignum (s) struct S_Bignum *s; {
  33.     Object big;
  34.     register size;
  35.  
  36.     big = Make_Uninitialized_Bignum (size = s->usize);
  37.     BIGNUM(big)->minusp = s->minusp;
  38.     BIGNUM(big)->usize = size;
  39.     bcopy ((char *)s->data, (char *)BIGNUM(big)->data,
  40.       size * sizeof (gran_t));
  41.     return big;
  42. }
  43.  
  44. Object Make_Bignum (buf, neg, radix) char *buf; {
  45.     Object big;
  46.     register char *p;
  47.     register c;
  48.     register size = (strlen (buf) + 4) / 4;
  49.  
  50.     big = Make_Uninitialized_Bignum (size);
  51.     BIGNUM(big)->minusp = neg ? True : False;
  52.     p = buf;
  53.     while (c = *p++) {
  54.     Bignum_Mult_In_Place (BIGNUM(big), radix);
  55.     if (radix == 16) {
  56.         if (isupper (c))
  57.         c = tolower (c);
  58.         if (c >= 'a')
  59.         c = '9' + c - 'a' + 1;
  60.     }
  61.     Bignum_Add_In_Place (BIGNUM(big), c - '0');
  62.     }
  63.     Bignum_Normalize_In_Place (BIGNUM(big)); /* to avoid -0 */
  64.     return big;
  65. }
  66.  
  67. Object Reduce_Bignum (x) Object x; {
  68.     register i;
  69.     register struct S_Bignum *p = BIGNUM(x);
  70.  
  71.     if (p->usize > 2 || (p->usize == 2 && p->data[1] >= 32768))
  72.     return x;
  73.     i = Bignum_To_Integer (x);
  74.     if (!FIXNUM_FITS(i))
  75.     return x;
  76.     return Make_Fixnum (i);
  77. }
  78.  
  79. Bignum_Mult_In_Place (x, n) register struct S_Bignum *x; {
  80.     register i = x->usize;
  81.     register gran_t *p = x->data;
  82.     register j;
  83.     register unsigned k = 0;
  84.     
  85.     for (j = 0; j < i; ++j) {
  86.     k += n * *p;
  87.         *p++ = k;
  88.     k >>= 16;
  89.     }
  90.     if (k) {
  91.     if (i >= x->size)
  92.         Panic ("Bignum_Mult_In_Place");
  93.     *p++ = k;
  94.     x->usize++;
  95.     }
  96. }
  97.  
  98. Bignum_Add_In_Place (x, n) register struct S_Bignum *x; {
  99.     register i = x->usize;
  100.     register gran_t *p = x->data;
  101.     register j = 0;
  102.     register unsigned k = n;
  103.  
  104.     if (i == 0) goto extend;
  105.     k += *p;
  106.     *p++ = k;
  107.     while (k >>= 16) {
  108.     if (++j >= i) {
  109.  extend:
  110.         if (i >= x->size)
  111.         Panic ("Bignum_Add_In_Place");
  112.         *p++ = k;
  113.         x->usize++;
  114.         return;
  115.     }
  116.     k += *p;
  117.     *p++ = k;
  118.     }
  119. }
  120.  
  121. Bignum_Div_In_Place (x, n) register struct S_Bignum *x; {
  122.     register i = x->usize;
  123.     register gran_t *p = x->data + i;
  124.     register unsigned k = 0;
  125.     for ( ; i; --i) {
  126.     k <<= 16;
  127.     k += *--p;
  128.     *p = k / n;
  129.     k %= n;
  130.     }
  131.     Bignum_Normalize_In_Place (x);
  132.     return k;
  133. }
  134.  
  135. Bignum_Normalize_In_Place (x) register struct S_Bignum *x; {
  136.     register i = x->usize;
  137.     register gran_t *p = x->data + i;
  138.     while (i && !*--p)
  139.     --i;
  140.     x->usize = i;
  141.     if (!i)
  142.     x->minusp = False;
  143. }
  144.  
  145. Print_Bignum (port, x) Object port, x; {
  146.     register char *p;
  147.     char *buf;
  148.     register size;
  149.     struct S_Bignum *big;
  150.     Alloca_Begin;
  151.     
  152.     if (Bignum_Zero (x)) {
  153.     Printf (port, "0");
  154.     return;
  155.     }
  156.     
  157.     size = BIGNUM(x)->usize * 5 + 3;
  158.     Alloca (buf, char*, size + 1);
  159.     p = buf + size;
  160.     *p = 0;
  161.  
  162.     size = (sizeof (struct S_Bignum) - sizeof (gran_t))
  163.     + BIGNUM(x)->usize * sizeof (gran_t);
  164.     Alloca (big, struct S_Bignum*, size);
  165.     bcopy ((char *)POINTER(x), (char *)big, size);
  166.     big->size = BIGNUM(x)->usize;
  167.  
  168.     while (big->usize) {
  169.     register unsigned bigdig = Bignum_Div_In_Place (big, 10000);
  170.     *--p = '0' + bigdig % 10;
  171.     bigdig /= 10;
  172.     *--p = '0' + bigdig % 10;
  173.     bigdig /= 10;
  174.     *--p = '0' + bigdig % 10;
  175.     bigdig /= 10;
  176.     *--p = '0' + bigdig;
  177.     }
  178.     while (*p == '0')
  179.     ++p;
  180.     if (Truep (BIGNUM(x)->minusp))
  181.     Printf (port, "-");
  182.     Format (port, p, strlen (p), 0, (Object *)0);
  183.     Alloca_End;
  184. }
  185.  
  186. Object Bignum_To_String (x, radix) Object x; {
  187.     register char *p;
  188.     char *buf;
  189.     register unsigned div, ndig, size;
  190.     struct S_Bignum *big;
  191.     Object ret;
  192.     Alloca_Begin;
  193.     
  194.     if (Bignum_Zero (x))
  195.     return Make_String ("0", 1);
  196.     
  197.     size = BIGNUM(x)->usize * (radix == 2 ? 17 : 6) + 3;
  198.     Alloca (buf, char*, size + 1);
  199.     p = buf + size;
  200.     *p = 0;
  201.  
  202.     size = (sizeof (struct S_Bignum) - sizeof (gran_t))
  203.     + BIGNUM(x)->usize * sizeof (gran_t);
  204.     Alloca (big, struct S_Bignum*, size);
  205.     bcopy ((char *)POINTER(x), (char *)big, size);
  206.     big->size = BIGNUM(x)->usize;
  207.  
  208.     switch (radix) {
  209.     case 2:
  210.     div = 65536; ndig = 16; break;
  211.     case 8:
  212.     div = 32768; ndig = 5; break;
  213.     case 10:
  214.     div = 10000; ndig = 4; break;
  215.     case 16:
  216.     div = 65536; ndig = 4; break;
  217.     }
  218.  
  219.     while (big->usize) {
  220.     register unsigned bigdig = Bignum_Div_In_Place (big, div);
  221.     register i;
  222.     for (i = 0; i < ndig; i++) {
  223.         *--p = '0' + bigdig % radix;
  224.         if (*p > '9')
  225.         *p = 'A' + (*p - '9') - 1;
  226.         bigdig /= radix;
  227.     }
  228.     }
  229.     while (*p == '0')
  230.     ++p;
  231.     if (Truep (BIGNUM(x)->minusp))
  232.     *--p = '-';
  233.     ret = Make_String (p, strlen (p));
  234.     Alloca_End;
  235.     return ret;
  236. }
  237.  
  238.  
  239. Bignum_To_Integer (x) Object x; {
  240.     unsigned n = 0;
  241.     int s = BIGNUM(x)->usize;
  242.  
  243.     if (s) {
  244.     n = BIGNUM(x)->data[0];
  245.     if (s > 1) {
  246.         n |= BIGNUM(x)->data[1] << 16;
  247.         if (s > 2)
  248. err:        
  249.         Primitive_Error ("integer out of range: ~s", x);
  250.     }
  251.     }
  252.     if (Truep (BIGNUM(x)->minusp)) {
  253.     if (n > (~(unsigned)0 >> 1) + 1)
  254.         goto err;
  255.     return -n;
  256.     } else {
  257.     if (n > ~(unsigned)0 >> 1)
  258.         goto err;
  259.     return n;
  260.     }
  261. }
  262.  
  263. Object Integer_To_Bignum (i) {
  264.     Object big = Make_Uninitialized_Bignum (2);
  265.     unsigned n = i;
  266.  
  267.     if (i < 0) {
  268.     BIGNUM(big)->minusp = True;
  269.     n = -i;
  270.     }
  271.     BIGNUM(big)->data[0] = n;
  272.     BIGNUM(big)->data[1] = n >> 16;
  273.     BIGNUM(big)->usize = 2;
  274.     Bignum_Normalize_In_Place (BIGNUM(big));
  275.     return big;
  276. }
  277.  
  278. Object Unsigned_To_Bignum (i) unsigned i; {
  279.     Object big = Make_Uninitialized_Bignum (2);
  280.  
  281.     BIGNUM(big)->data[0] = i;
  282.     BIGNUM(big)->data[1] = i >> 16;
  283.     BIGNUM(big)->usize = 2;
  284.     Bignum_Normalize_In_Place (BIGNUM(big));
  285.     return big;
  286. }
  287.  
  288. Object Double_To_Bignum (d) double d; {         /* Truncates the double */
  289.     Object big;
  290.     int expo, size;
  291.     double mantissa = frexp (d, &expo);
  292.     register gran_t *p;
  293.     
  294.     if (expo <= 0 || mantissa == 0.0)
  295.     return Make_Uninitialized_Bignum (0);
  296.     size = (expo + (16-1)) / 16;
  297.     big = Make_Uninitialized_Bignum (size);
  298.     BIGNUM(big)->usize = size;
  299.     if (mantissa < 0.0) {
  300.     BIGNUM(big)->minusp = True;
  301.     mantissa = -mantissa;
  302.     }
  303.     p = BIGNUM(big)->data;
  304.     bzero ((char *)p, size * sizeof (gran_t));
  305.     p += size;
  306.     if (expo &= (16-1))
  307.     mantissa = ldexp (mantissa, expo - 16);
  308.     while (mantissa != 0.0) {
  309.     if (--size < 0)
  310.         break;        /* inexact */
  311.     mantissa *= 65536.0;
  312.     *--p = (int)mantissa;
  313.     mantissa -= *p;
  314.     }
  315.     Bignum_Normalize_In_Place (BIGNUM(big)); /* Probably not needed */
  316.     return Reduce_Bignum (big);
  317. }
  318.  
  319. double Bignum_To_Double (x) Object x; {   /* error if it ain't fit */
  320.     double rx = 0.0;
  321.     register i = BIGNUM(x)->usize;
  322.     register gran_t *p = BIGNUM(x)->data + i;
  323.  
  324.     for (i = BIGNUM(x)->usize; --i >= 0; ) {
  325.     if (rx >= HUGE / 65536.0)
  326.         Primitive_Error ("cannot coerce to real: ~s", x);
  327.     rx *= 65536.0;
  328.     rx += *--p;
  329.     }
  330.     if (Truep (BIGNUM(x)->minusp))
  331.     rx = -rx;
  332.     return rx;
  333. }
  334.  
  335. Bignum_Zero (x) Object x; {
  336.     return BIGNUM(x)->usize == 0;
  337. }
  338.  
  339. Bignum_Negative (x) Object x; {
  340.     return Truep (BIGNUM(x)->minusp);
  341. }
  342.  
  343. Bignum_Positive (x) Object x; {
  344.     return !Truep (BIGNUM(x)->minusp) && BIGNUM(x)->usize != 0;
  345. }
  346.  
  347. Bignum_Even (x) Object x; {
  348.     return BIGNUM(x)->usize == 0 || (BIGNUM(x)->data[0] & 1) == 0;
  349. }
  350.  
  351. Object Bignum_Abs (x) Object x; {
  352.     Object big;
  353.  
  354.     big = Copy_Bignum (x);
  355.     BIGNUM(big)->minusp = False;
  356.     return big;
  357. }
  358.  
  359. Bignum_Mantissa_Cmp (x, y) register struct S_Bignum *x, *y; {
  360.     register i = x->usize;
  361.     if (i < y->usize)
  362.     return -1;
  363.     else if (i > y->usize)
  364.     return 1;
  365.     else {
  366.     register gran_t *xbuf = x->data + i;
  367.     register gran_t *ybuf = y->data + i;
  368.     for ( ; i; --i) {
  369.         register n;
  370.         if (n = (int)*--xbuf - (int)*--ybuf)
  371.         return n;
  372.     }
  373.     return 0;
  374.     }
  375. }
  376.  
  377. Bignum_Cmp (x, y) register struct S_Bignum *x, *y; {
  378.     register xm = Truep (x->minusp);
  379.     register ym = Truep (y->minusp);
  380.     if (xm) {
  381.     if (ym)
  382.         return -Bignum_Mantissa_Cmp (x, y);
  383.     else return -1;
  384.     } else {
  385.     if (ym)
  386.         return 1;
  387.     else return Bignum_Mantissa_Cmp (x, y);
  388.     }
  389. }
  390.  
  391. Bignum_Equal (x, y) Object x, y; {
  392.     return Bignum_Cmp (BIGNUM(x), BIGNUM(y)) == 0;
  393. }
  394.  
  395. Bignum_Less (x, y) Object x, y; {
  396.     return Bignum_Cmp (BIGNUM(x), BIGNUM(y)) < 0;
  397. }
  398.  
  399. Bignum_Greater (x, y) Object x, y; {
  400.     return Bignum_Cmp (BIGNUM(x), BIGNUM(y)) > 0;
  401. }
  402.  
  403. Bignum_Eq_Less (x, y) Object x, y; {
  404.     return Bignum_Cmp (BIGNUM(x), BIGNUM(y)) <= 0;
  405. }
  406.  
  407. Bignum_Eq_Greater (x, y) Object x, y; {
  408.     return Bignum_Cmp (BIGNUM(x), BIGNUM(y)) >= 0;
  409. }
  410.  
  411. Object General_Bignum_Plus_Minus (x, y, neg) Object x, y; {
  412.     Object big;
  413.     int size, xsize, ysize, xminusp, yminusp;
  414.     GC_Node2;
  415.  
  416.     GC_Link2 (x,y);
  417.     xsize = BIGNUM(x)->usize;
  418.     ysize = BIGNUM(y)->usize;
  419.     xminusp = Truep (BIGNUM(x)->minusp);
  420.     yminusp = Truep (BIGNUM(y)->minusp);
  421.     if (neg)
  422.     yminusp = !yminusp;
  423.     size = xsize > ysize ? xsize : ysize;
  424.     if (xminusp == yminusp)
  425.     size++;
  426.     big = Make_Uninitialized_Bignum (size);
  427.     BIGNUM(big)->usize = size;
  428.     GC_Unlink;
  429.  
  430.     if (xminusp == yminusp) {
  431.     /* Add x and y */
  432.     register unsigned k = 0;
  433.     register i;
  434.     register gran_t *xbuf = BIGNUM(x)->data;
  435.     register gran_t *ybuf = BIGNUM(y)->data;
  436.     register gran_t *zbuf = BIGNUM(big)->data;
  437.     for (i = 0; i < size; ++i) {
  438.         if (i < xsize)
  439.         k += *xbuf++;
  440.         if (i < ysize)
  441.         k += *ybuf++;
  442.         *zbuf++ = k;
  443.         k >>= 16;
  444.     }
  445.     } else {
  446.     if (Bignum_Mantissa_Cmp (BIGNUM(x), BIGNUM(y)) < 0) {
  447.         Object temp = x;
  448.         x = y; y = temp;
  449.         xsize = ysize;
  450.         ysize = BIGNUM(y)->usize;
  451.         xminusp = yminusp;
  452.     }
  453.     /* Subtract y from x */
  454.     {
  455.     register unsigned k = 1;
  456.     register i;
  457.     register gran_t *xbuf = BIGNUM(x)->data;
  458.     register gran_t *ybuf = BIGNUM(y)->data;
  459.     register gran_t *zbuf = BIGNUM(big)->data;
  460.     for (i = 0; i < size; ++i) {
  461.         if (i < xsize)
  462.         k += *xbuf++;
  463.         else Panic ("General_Bignum_Plus_Minus");
  464.         if (i < ysize)
  465.         k += ~*ybuf++ & 0xFFFF;
  466.         else k += 0xFFFF;
  467.         *zbuf++ = k;
  468.         k >>= 16;
  469.     }
  470.     }
  471.     }
  472.     BIGNUM(big)->minusp = xminusp ? True : False;
  473.     Bignum_Normalize_In_Place (BIGNUM(big));
  474.     return Reduce_Bignum (big);
  475. }
  476.  
  477. Object Bignum_Plus (x, y) Object x, y; {   /* bignum + bignum */
  478.     return General_Bignum_Plus_Minus (x, y, 0);
  479. }
  480.  
  481. Object Bignum_Minus (x, y) Object x, y; {   /* bignum - bignum */
  482.     return General_Bignum_Plus_Minus (x, y, 1);
  483. }
  484.  
  485. Object Bignum_Fixnum_Multiply (x, y) Object x, y; {   /* bignum * fixnum */
  486.     Object big;
  487.     register size, xsize, i;
  488.     register gran_t *xbuf, *zbuf;
  489.     int yn = FIXNUM(y);
  490.     register unsigned yl, yh;
  491.     GC_Node;
  492.     
  493.     GC_Link (x);
  494.     xsize = BIGNUM(x)->usize;
  495.     size = xsize + 2;
  496.     big = Make_Uninitialized_Bignum (size);
  497.     BIGNUM(big)->usize = size;
  498.     if (Truep (BIGNUM(x)->minusp) != (yn < 0))
  499.     BIGNUM(big)->minusp = True;
  500.     bzero ((char *)BIGNUM(big)->data, size * sizeof (gran_t));
  501.     xbuf = BIGNUM(x)->data;
  502.     if (yn < 0)
  503.     yn = -yn;
  504.     yl = yn & 0xFFFF;
  505.     yh = yn >> 16;
  506.     zbuf = BIGNUM(big)->data;
  507.     for (i = 0; i < xsize; ++i) {
  508.     register unsigned xf = xbuf[i];
  509.     register unsigned k = 0;
  510.     register gran_t *r = zbuf + i;
  511.     k += xf * yl + *r;
  512.     *r++ = k;
  513.     k >>= 16;
  514.     k += xf * yh + *r;
  515.     *r++ = k;
  516.     k >>= 16;
  517.     *r = k;
  518.     }
  519.     GC_Unlink;
  520.     Bignum_Normalize_In_Place (BIGNUM(big));
  521.     return Reduce_Bignum (big);
  522. }
  523.  
  524. Object Bignum_Multiply (x, y) Object x, y; {   /* bignum * bignum */
  525.     Object big;
  526.     register size, xsize, ysize, i, j;
  527.     register gran_t *xbuf, *ybuf, *zbuf;
  528.     GC_Node2;
  529.     
  530.     GC_Link2 (x, y);
  531.     xsize = BIGNUM(x)->usize;
  532.     ysize = BIGNUM(y)->usize;
  533.     size = xsize + ysize;
  534.     big = Make_Uninitialized_Bignum (size);
  535.     BIGNUM(big)->usize = size;
  536.     if (!EQ(BIGNUM(x)->minusp, BIGNUM(y)->minusp))
  537.     BIGNUM(big)->minusp = True;
  538.     bzero ((char *)BIGNUM(big)->data, size * sizeof (gran_t));
  539.     xbuf = BIGNUM(x)->data;
  540.     ybuf = BIGNUM(y)->data;
  541.     zbuf = BIGNUM(big)->data;
  542.     for (i = 0; i < xsize; ++i) {
  543.     register unsigned xf = xbuf[i];
  544.     register unsigned k = 0;
  545.     register gran_t *p = ybuf;
  546.     register gran_t *r = zbuf + i;
  547.     for (j = 0; j < ysize; ++j) {
  548.         k += xf * *p++ + *r;
  549.         *r++ = k;
  550.         k >>= 16;
  551.     }
  552.     *r = k;
  553.     }
  554.     GC_Unlink;
  555.     Bignum_Normalize_In_Place (BIGNUM(big));
  556.     return Reduce_Bignum (big);
  557. }
  558.  
  559. /* Returns cons cell (quotient . remainder); cdr is a fixnum
  560.  */
  561. Object Bignum_Fixnum_Divide (x, y) Object x, y; {   /* bignum / fixnum */
  562.     Object big;
  563.     register xsize, i;
  564.     register gran_t *xbuf, *zbuf;
  565.     int yn = FIXNUM(y);
  566.     int xminusp, yminusp = 0;
  567.     register unsigned rem;
  568.     GC_Node;
  569.  
  570.     GC_Link (x);
  571.     if (yn < 0) {
  572.     yn = -yn;
  573.     yminusp = 1;
  574.     }
  575.     if (yn > 0xFFFF) {
  576.     big = Integer_To_Bignum (FIXNUM(y));
  577.     GC_Unlink;
  578.     return Bignum_Divide (x, big);
  579.     }
  580.     xsize = BIGNUM(x)->usize;
  581.     big = Make_Uninitialized_Bignum (xsize);
  582.     BIGNUM(big)->usize = xsize;
  583.     xminusp = Truep (BIGNUM(x)->minusp);
  584.     if (xminusp != yminusp)
  585.     BIGNUM(big)->minusp = True;
  586.     xbuf = BIGNUM(x)->data;
  587.     zbuf = BIGNUM(big)->data;
  588.     rem = 0;
  589.     for (i = xsize; --i >= 0; ) {
  590.     rem <<= 16;
  591.     rem += xbuf[i];
  592.     zbuf[i] = rem / yn;
  593.     rem %= yn;
  594.     }
  595.     GC_Unlink;
  596.     Bignum_Normalize_In_Place (BIGNUM(big));
  597.     if (xminusp)
  598.     rem = -(int)rem;
  599.     return Cons (Reduce_Bignum (big), Make_Fixnum ((int)rem));
  600. }
  601.  
  602. /* Returns cons cell (quotient . remainder); cdr is a fixnum
  603.  */
  604. Object Bignum_Divide (x, y) Object x, y; {   /* bignum / bignum */
  605.     struct S_Bignum *dend, *dor;
  606.     int quotsize, dendsize, dorsize, scale;
  607.     unsigned dor1, dor2;
  608.     Object quot, rem;
  609.     register gran_t *qp, *dendp;
  610.     GC_Node2;
  611.     Alloca_Begin;
  612.  
  613.     if (BIGNUM(y)->usize < 2)
  614.     return Bignum_Fixnum_Divide (x, Make_Fixnum (Bignum_To_Integer (y)));
  615.  
  616.     GC_Link2 (x, y);
  617.     quotsize = BIGNUM(x)->usize - BIGNUM(y)->usize + 1;
  618.     if (quotsize < 0)
  619.     quotsize = 0;
  620.     quot = Make_Uninitialized_Bignum (quotsize);
  621.     GC_Unlink;
  622.  
  623.     dendsize = (sizeof (struct S_Bignum) - sizeof (gran_t))
  624.     + (BIGNUM(x)->usize + 1) * sizeof (gran_t);
  625.     Alloca (dend, struct S_Bignum*, dendsize);
  626.     bcopy ((char *)POINTER(x), (char *)dend, dendsize);
  627.     dend->size = BIGNUM(x)->usize + 1;
  628.  
  629.     if (quotsize == 0 || Bignum_Mantissa_Cmp (dend, BIGNUM(y)) < 0)
  630.     goto zero;
  631.  
  632.     dorsize = (sizeof (struct S_Bignum) - sizeof (gran_t))
  633.     + BIGNUM (y)->usize * sizeof (gran_t);
  634.     Alloca (dor, struct S_Bignum*, dorsize);
  635.     bcopy ((char *)POINTER(y), (char *)dor, dorsize);
  636.     dor->size = dorsize = BIGNUM(y)->usize;
  637.  
  638.     scale = 65536 / (dor->data[dor->usize - 1] + 1);
  639.     Bignum_Mult_In_Place (dend, scale);
  640.     if (dend->usize < dend->size)
  641.     dend->data[dend->usize++] = 0;
  642.     Bignum_Mult_In_Place (dor, scale);
  643.  
  644.     BIGNUM(quot)->usize = BIGNUM(quot)->size;
  645.     qp = BIGNUM(quot)->data + BIGNUM(quot)->size;
  646.     dendp = dend->data + dend->usize;
  647.     dor1 = dor->data[dor->usize - 1];
  648.     dor2 = dor->data[dor->usize - 2];
  649.         
  650.     while (qp > BIGNUM(quot)->data) {
  651.     unsigned msw, guess;
  652.     int k;
  653.     register gran_t *dep, *dop, *edop;
  654.     
  655.     msw = dendp[-1] << 16 | dendp[-2];
  656.     guess = msw / dor1;
  657.     if (guess >= 65536)    /* [65535, 0, 0] / [65535, 65535] */
  658.         guess = 65535;
  659.     for (;;) {
  660.         unsigned d1, d2, d3;
  661.         d3 = dor2 * guess;
  662.         d2 = dor1 * guess + (d3 >> 16);
  663.         d3 &= 0xFFFF;
  664.         d1 = d2 >> 16;
  665.         d2 &= 0xFFFF;
  666.         if (d1 < dendp[-1] || (d1 == dendp[-1] &&
  667.                    (d2 < dendp[-2] || (d2 == dendp[-2] &&
  668.                                d3 <= dendp[-3]))))
  669.         break;
  670.         --guess;
  671.     }
  672.     --dendp;
  673.     k = 0;
  674.     dep = dendp - dorsize;
  675.     for (dop = dor->data, edop = dop + dor->usize; dop < edop; ) {
  676.         register unsigned prod = *dop++ * guess;
  677.         k += *dep;
  678.         k -= prod & 0xFFFF;
  679.         *dep++ = k;
  680.         k >>= 16;
  681.         k -= prod >> 16;
  682.     }
  683.     k += *dep;
  684.     *dep = k;
  685.     if (k < 0) {
  686.         k = 0;
  687.         dep = dendp - dorsize;
  688.         for (dop = dor->data, edop = dop + dor->usize; dop < edop; ) {
  689.         k += *dep + *dop++;
  690.         *dep++ = k;
  691.         k >>= 16;
  692.         }
  693.         k += *dep;
  694.         *dep = k;
  695.         --guess;
  696.     }
  697.     *--qp = guess;
  698.     }
  699.     
  700.     if (Bignum_Div_In_Place (dend, scale))
  701.     Panic ("Bignum_Div scale");
  702.  zero:
  703.     if (Truep (dend->minusp = BIGNUM(x)->minusp) != Truep (BIGNUM(y)->minusp))
  704.     BIGNUM(quot)->minusp = True;
  705.     Bignum_Normalize_In_Place (BIGNUM(quot));
  706.     Bignum_Normalize_In_Place (dend);
  707.     GC_Link (quot);
  708.     rem = Reduce_Bignum (Copy_S_Bignum (dend));
  709.     GC_Unlink;
  710.     Alloca_End;
  711.     return Cons (Reduce_Bignum (quot), rem);
  712. }
  713.